### PPP
### O-Logits

rm(list = ls())

### libraries and packages
pkg <- c("plyr", "dplyr", "tidyr", 
         "MASS", "multiwayvcov", 
         "stargazer", "lmtest", "doBy", "ggplot2", 
         "mediation", "vars", "DataCombine", "foreign",
         "mgcv", "ordinal", "reshape2", 
         "xtable", "sandwich", "rms")
lapply(pkg, require, character.only = TRUE)

# rid of E
options(scipen=999)
### setting wd for inputs
setwd("~/Dropbox/Pollution and Public Perceptions in China/Data and Analysis/Data")

# loading the data
load("ppp_cleaned")

#####################################
#####################################
## o-logits
#####################################
#####################################

l_sato <- lrm(factor(l_sat) ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ + kids + married + income + gender + ccp + parade + national + age, p2, y=TRUE, x = TRUE)
robcov(l_sato, cluster=p2$day) 
coeftest(l_sato, robcov(l_sato, cluster = p2$day)$var)
seo1 <- coeftest(l_sato, robcov(l_sato, cluster = p2$day)$var)

colnames(p2)
c_sato <- lrm(factor(c_sat) ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ  + kids + married + income + gender + ccp + parade + national + age, p2, y=TRUE, x = TRUE)
robcov(c_sato, cluster=p2$day) 
coeftest(c_sato, robcov(c_sato, cluster = p2$day)$var)
seo2 <- coeftest(c_sato, robcov(c_sato, cluster = p2$day)$var)

l_watcho <- lrm(factor(l_watch) ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ + kids + married + income + gender + ccp + parade + national + age, p2, y=TRUE,x = TRUE)
robcov(l_watcho, cluster=p2$day)
coeftest(l_watcho, robcov(l_watcho, cluster = p2$day)$var)
seo3 <- coeftest(l_watcho, robcov(l_watcho, cluster = p2$day)$var)

c_watcho <- lrm(factor(c_watch) ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ  + kids + married + income + gender + ccp + parade + national + age, p2, y=TRUE, x = TRUE)
robcov(c_watcho, cluster=p2$day) 
coeftest(c_watcho, robcov(c_watcho, cluster = p2$day)$var)
seo4 <- coeftest(c_watcho, robcov(c_watcho, cluster = p2$day)$var)

anti_westo <- lrm(factor(anti_west) ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ + kids + married + income + gender + ccp + parade + national + age, p2, y=TRUE,x = TRUE)
robcov(anti_westo, cluster=p2$day) # robust to clustered standard errors
coeftest(anti_westo, robcov(anti_westo, cluster = p2$day)$var)
seo6 <- coeftest(anti_westo, robcov(anti_westo, cluster = p2$day)$var)

day_polluteo <- lrm(factor(day_pollute) ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ + kids + married + income + gender + ccp + parade + national + age, p2, y=TRUE,x = TRUE)
robcov(day_polluteo, cluster=p2$day) 
coeftest(day_polluteo, robcov(day_polluteo, cluster = p2$day)$var)
seo7 <- coeftest(day_polluteo, robcov(day_polluteo, cluster = p2$day)$var)

life_sato <- lrm(factor(life_sat) ~ aqi + aqi.lagged + aqi.lagged2 + hukou + insider + soe + educ + kids + married + income + gender + ccp + parade + national + age, p2, y=TRUE,x = TRUE)
robcov(life_sato, cluster=p2$day) 
coeftest(life_sato, robcov(life_sato, cluster = p2$day)$var)
seo8 <- coeftest(life_sato, robcov(life_sato, cluster = p2$day)$var)

# Making ologit table

fit.list<- list(l_sato, c_sato, l_watcho, c_watcho, anti_westo)
se.list <- list(seo1[,2], seo2[,2], seo3[, 2], seo4[, 2], seo6[, 2])
column.labels <- c("Local Government", "Central Government", "Local Oversight", "Central Oversight","Anti-West")
cov.labs.c <-c("AQI", "AQI Lagged", "AQI Lagged 2","Hukou Status", "Regime Insider","SOE Employee", "Education","Children","Married","Income", "Female", "CCP Member","Parade Period","Holiday","Age")

## TABLE A5
sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/ologit.tex")
stargazer(fit.list,
          se = se.list,
          font.size = "large", digits = 4, 
          column.labels = column.labels, 
          covariate.labels = cov.labs.c,
          dep.var.labels.include = FALSE,
          model.numbers = FALSE,
          style = "qje", no.space = FALSE, 
          keep.stat = c("n"),
          star.cutoffs = c(.10, .05, .01), 
          omit = "y",
          header = F,
          notes = "robust standard errors clustered at the day level in parentheses", 
          float = F)
sink()
